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target  identification  and  weapon  guidance. 

Practical  implementation  of  feature  extracting 
algorithms  is  constrained  by  size,  weight,  power  and  data 
throughput  requirements  imposed  by  the  particular 
application.  These  constraints  limit  the  amount  of 
computation  that  can  be  accomplished  on  the  input  image 
data.  In  order  to  allow  more  of  the  available  computational 
time  and  hardware  to  be  devoted  to  sophisticated  and 
computationally  intensive  feature  classification  routines, 
it  is  desirable  to  minimize  the  preliminary  feature 
extracting  calculations. 

One  approach  to  the  problem  of  reducing  calculations  is 
to  map  the  original  intensity  image  into  a  smaller  image 
before  applying  feature  extracting  algorithms.  Several  such 
image  data  reducing  techniques  were  implemented  in  [1], 
followed  by  standard  (Roberts  [2-4]  and  Kirsch  [2,3])  edge 
detecting  operations  on  the  reduced  images.  Since  edges  in 
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an  Image  contain  much  of  the  information  necessary  to 
classify  objects,  edge  detectability  and  edge  quality  were 
used  in  that  work  as  measures  of  feature  information  loss. 
The  outputs  of  the  standard  edge  detectors  were 
quantitatively  evaluated  to  determine  the  effect  of  the 
image  reducing  techniques  on  edge  content. 

For  the  image  reducing  techniques  implemented  in  [1], 
it  is  evident  that  the  straightforward  application  of 
standard  edge  detectors  to  the  reduced  images  does  not  fully 
extract  the  edge  information  that  is  available.  This  is 
demonstrated  by  the  substantially  better  results  achieved  by 
applying  the  Kirsch  operator  to  re-expanded  versions  of  the 
reduced  images.  The  re-expansion  was  accomplished  by  simply 
duplicating  reduced  image  picture  elements  (pixels),  so  the 
information  content  would  seem  to  be  unaffected.  The 
two-step  sequence  of  re-expansion  followed  by  Kirsch  edge 
detection  is  equivalent  to  a  new  edge  detection  scheme 
operating  on  a  reduced  image  and  generating  multiple  output 
pixels  for  each  input  pixel.  The  re-expansion  operation 
Inherently  Incorporates  knowledge  of  the  original  image 
size,  adding  information  that  was  not  available  in  the 
reduced  image.  However,  the  two-step  sequence  does  not 
incorporate  any  knowledge  of  the  reducing  technique  used. 

The  purpose  of  this  research  Is  to  find  techniques  for 
extracting  edges  from  reduced  images  based  on  knowledge  of 
the  specific  image  reducing  techniques  used.  Incorporating 
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knowledge  of  the  Image  reducing  techniques  used  should 
result  in  edge  outputs  that  are  more  representative  of  the 
edge  information  available  in  reduced  images. 

The  search  for  better  ways  to  extract  edges  from 
reduced  images  leads  naturally  to  consideration  of  how  the 
image  data  should  be  reduced  in  the  first  place.  As  a 
result  this  investigation  centers  on  finding  methods  for 
reducing  an  image  so  as  to  maximize  the  retention  of  edge 
information  which  is  subsequently  extracted.  Since  lower 
mean-square  error  (MSE)  in  the  intensity  domain  was  usually 
(with  one  exception)  found  to  correspond  to  lower  MSE  in  the 
edge  domain  in  [1],  The  Hotelling  transform  [2,5]  (also 
called  the  discrete  Karhunen-Loeve  transform  [2,4])  is 
initially  investigated.  The  Hotelling  transform  minimizes 
the  intensity  mean-square  error  (IMSE)  between  an  original 
image  and  one  that  has  been  reconstructed  by  an  inverse 
Hotelling  transformation  using  partial  transform 
coefficients  (i.e.  reduced  image  data).  Subsequently  a 
measure  of  edge  loss  called  the  gradient  mean-square  error 
(GMSE)  based  on  the  Roberts  gradient  is  defined,  and  the 
relationship  between  GMSE  and  IMSE  is  derived.  Finally,  a 
linear  transformation  is  introduced  which  reduces  an  image 
based  on  minimizing  GMSE  (and  in  that  sense  maximizing  edge 
retention) . 
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2.  OBJECTIVE 


The  primary  purpose  of  this  research  is  to  find 
techniques  for  extracting  edges  from  reduced  images  based  on 
knowledge  of  the  specific  image  data  reducing  techniques 
sed . 
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3.  EDGE  EVALUATION  TECHNIQUES 


Several  techniques  are  used  In  subsequent  sections  of 
this  report  to  compare  edges  In  original  and  reconstructed 
Images.  Three  original  256  by  256  Images  differing  In  edge 
content  are  used  In  this  study.  They  are  shown  In  Figure  I 
and  will  be  referred  to  as  Image  A  (upper  left).  Image  B 
(upper  right),  and  Image  C  (lower  left). 

It  is  difficult  to  find  edge  evaluation  techniques 
applicable  to  real-world  imagery.  This  is  due  in  part  to 
the  inability  to  clearly  define  the  actual  location  of  edges 
in  an  image  and  thus  form  a  reliable  basis  for  comparison. 
However,  a  paper  by  Kitchen  and  Rosenfeld  [4]  proposes  an 
edge  quality  evaluation  based  on  edge  coherence 
incorporating  connectedness,  thinness,  and  gradient 
directions  as  seen  by  3x3  neighborhoods.  The  specific 
implementation  used  here  is  described  in  [I].  In  this 
implementation  a  Kirsch  edge  detector  is  applied  to  an 
image.  Each  pixel  location  in  the  resulting  Kirsch  image 
that  has  a  Kirsch  response  equal  to  or  greater  than  a 
selected  threshold  is  assigned  an  edge  coherence  value 

E  -  wC  +  (l-w)T  <3-D 
where  C  and  T  are  measures  of  edge  continuity  and  thinness, 
respectively,  and  w  is  a  weighting  factor.  E,  C,  T,  &  w  are 
constrained  to  lie  in  the  interval  [0,1].  In  this  study  w  - 
0.8  as  recommended  in  [6].  The  average  value  of  E  for  a 
selected  threshold  is  used  as  an  overall  measure  of  image 
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edge  quality. 

In  order  to  determine  the  edge  retaining  abilities  of 
an  image  data  reducing  technique  some  measure  of  the  edge 
information  lost  in  the  transformation  is  needed.  One  such 
measure  based  on  the  Roberts  edge  detector  [2-5]  is  obtained 
by  calculating  a  gradient  vector  G  at  each  pixel  in  the 
original  and  reconstructed  images  and  determining  the 
mean-square  error  between  the  original  and  reconstructed 
gradient  vectors.  This  gradient  mean-square  error  (GMSE)  is 
used  as  a  measure  of  edge  information  loss.  For  a  pixel 
X(p,q)  at  image  location  ( r ow , co lumn ) =( p , q )  the  specific 
gradient  used  was  defined  as 

G=0.5*I*[(X(p,q-l)  +  X( p , q ) )  -  (X(p-l,q-l)  +  X(p-l,q))] 
+0. 5*  J* [ (X(p-l,q)  +  X( p , q ) )  -  (X(p-l,q-l)  +  X(p,q-1))] 

(3.2) 

where  I  and  J  are  unit  vectors  in  the  directions  down  and 
right,  respectively,  and  *  represents  multiplication. 

Two  other  MSE  based  measures  were  obtained  by  operating 
on  original  and  reconstructed  images  using  the  Roberts  (or 
Kirsch)  edge  detector  and  calculating  the  resulting  MSE 
between  the  original  and  reconstructed  edge  images  averaged 
over  all  pixel  locations  to  obtain  the  Roberts  (or  Kirsch) 
MSE.  A  fourth  MSE  measure  which  will  be  referred  to  as  the 


thresholded  Kirsch  MSE  is  obtained  by  calculating  the  Kirsch 
MSE  averaged  over  only  those  pixel  locations  where  the 


Kirsch  response  is  equal  to  or  greater  than  a  selected 
threshold . 

The  last  edge  comparison  technique  implemented  uses  an 
approach  from  signal  theory  which  calculates  the 
probabilities  of  edge  detection  PD  and  false  alarm  PF . 
These  calculations  use  binary  reference  and  test  images 
generated  by  globally  thresholding  the  edge  detector  output 
images.  The  resulting  edge  pixels  are  compared  to  give  PD 
and  PF  where: 

PD  is  the  fraction  of  edge  pixels  in  the  reference  image 
that  are  correctly  classified  as  edge  pixels  in  the  test 
image ; 

and 

PF  is  the  fraction  of  non-edge  pixels  in  the  reference  image 
that  are  incorrectly  classified  as  edge  pixels  in  the  test 
image . 
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4. 


IMAGE  DATA  REDUCTION  USING  THE  HOTELLING  TRANSFORM 


Since  lower  MSE  in  the  intensity  domain  was  found  in 
[1]  to  strongly  correlate  with  lower  MSE  in  the  edge  domain, 
the  effect  of  the  Hotelling  transform  on  edge  retention  is 


studied 

here  . 

The  Hotelling  transform 

minimi zes 

the 

IMSE 

between 

an 

original  image  and 

one  that 

has 

been 

reconstructed 

by  an  inverse  Hotelling 

t  r ans  f  o  rma  t i on 

using 

partial 

transform  coefficients. 

As  in 

[1]  , 

this 

investigation  concentrates  on  4:1  data  compression. 

The  Hotelling  transform  was  implemented  by  partitioning 
an  image  into  adjacent  non-overlapping  N  by  N  blocks.  Each 
N  by  N  block  is  row  scanned  to  form  a  vector  X  with  N**2 
elements,  where  **  represents  exponentiation.  Using 
notation  similar  to  that  in  [5],  a  complete  set  of  N**2 
Hotelling  transform  coefficients  forming  a  vector  Y  could  be 
calculated  using 

Y  =  A ( X  -  MX)  (4.1) 
where  the  vector  MX  is  the  ensemble  mean  of  the  vector  X 
averaged  over  all  blocks  in  the  image,  and  A  is  an  N**2  by 
N**2  matrix  whose  rows  are  the  eigenvectors  of  the  image 
covariance  matrix 

CX  =  E  {(X  -  MX)  (X  -  MX)'}  (4.2) 
where  E  is  the  expectation  operator  averaged  over  all  blocks 
in  the  image  and  the  prime  (')  indicates  transposition.  To 
achieve  image  data  reduction  ISAVE  coefficients  are 
generated  (ISAVE  <  N**2)  using  an  ISAVE  by  N**2  matrix  A. 


For  a  specified  amount  of  data  compression  the  Hotelling 
transform  minimizes  the  IMSE  in  the  reconstructed  image  by 
choosing  the  rows  of  A  to  be  the  ISAVE  eigenvectors 
corresponding  to  the  ISAVE  largest  eigenvalues  of  the 
covariance  matrix  CX.  To  achieve  4:1  image  data  reduction 
ISAVE  =  (N**2)/4  was  used.  For  each  N  by  N  block  a 
transform  coefficient  vector  Y  with  dimension  ISAVE  is 
generated.  Images  are  reconstructed  a  block  at  a  time  from 
the  reduced  image  data  using  the  inverse  Hotelling 
t  rans  f  orma t ion 

XE  ST  -  A '  Y  +  MX  (4.3) 
where  XEST  is  the  reconstructed  estimate  of  X. 

For  various  values  of  N  the  test  images  were  data 
compressed  and  reconstructed  using  the  Hotelling  transform. 
Then  the  evaluation  techniques  described  in  the  previous 
section  were  applied  to  the  reconstructed  images. 


For 

comparison 

purposes 

these 

same 

evaluat ion 

techniques 

were  also 

applied  to 

images 

which 

had 

been 

reduced 

by  simply 

averaging  2 

by  2 

blocks 

and 

then 

reconstructed  by  pixel  duplication.  The  results  are  shown 
in  Tables  1  &  2.  Table  1  values  were  calculated  for  all 
three  test  images.  Table  2  values  were  calculated  for  Image 
A  only  using  a  threshold  value  of  99  which  corresponds  to  a 
maximum  average  value  of  edge  coherence  E  for  the  original 
Image  A.  The  column  labelled  CDE  in  Table  1  will  be  defined 
in  the  next  section. 
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As  is  well  known,  for  any  specified  amount  of 
compression  the  reconstructed  IMSE  using  the  Hotelling 
transform  tends  to  decrease  as  the  block  size  increases. 
This  trend  is  evident  in  Table  1.  The  Hotelling  transform 
accomplishes  this  by  taking  advantage  of  correlations  that 
exist  between  neighboring  pixels.  As  N  gets  larger  than  the 
maximum  distance  between  pixels  with  significant 
correlation,  a  plot  of  IMSE  versus  N  would  tend  to  level 
off.  From  Table  1  there  appears  to  still  be  significant 
correlation  between  pixels  separated  by  a  distance  of  8 
since  IMSE  is  still  dropping  rapidly  for  N  =  8.  The  simple 
averaging  of  2  by  2  blocks  is  unable  to  take  advantage  of 
correlations  of  pixels  separated  by  distances  greater  than 
one.  Therefore,  as  expected,  the  Hotelling  transform  offers 
substantially  better  IMSE  performance  than  simple  averaging. 
Notably,  the  edge  evaluating  measures  shown  in  Tables  1  &  2 
also  demonstrate  that  the  Hotelling  transform  has 
significantly  better  edge  retaining  ability  than  simple 
averaging. 

A  major  obstacle  to  the  use  of  the  Hotelling  transform 
in  real-time  encoding  environments  is  the  calculation  of  the 
eigenvectors  of  the  block  covariance  matrix  CX  for  each 
image.  It  may  however  be  possible  to  generate  a  single 
Hotelling  transform  matrix  A  based  on  some  selected 
covariance  matrix  that  will  retain  edges  satisfactorily  over 
a  wide  variety  of  images.  If  necessary  a  small  set  of 


10 


selectable  A  matrices  might  be  stored  for  real-time  use.  To 
study  the  sensitivity  of  edge  retention  to  the  particular 
transformation  matrix  A  used,  the  test  Images  were  reduced 
and  reconstructed  using  Hotelling  transformations  based  on 
each  others  covariance  matrices.  The  resulting  GMSE  values 
are  given  in  Table  3.  Images  B  and  C  appear  quite 
insensitive  to  the  specific  covariance  matrix  used.  Image 
A,  which  contains  the  most  edges.  Is  the  most  sensitive  to 
the  covariance  matrix  selected.  Covariance  matrix  selection 
based  on  the  type  of  terrain  being  viewed  should  be 
practical . 
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5.  RELATIONSHIP  BETWEEN  GMSE  &  IMSE 

As  noted  in  the  previous  section,  there  appears  to  be  a 
strong  correlation  between  IMSE  and  the  edge  domain 
measures.  In  this  section  the  edge  loss  measure  GMSE  is 
found  to  be  related  to  the  IMSE  by 

GMSE  =  2*  ( IMSE  -  CDE )  (5.1) 

where  CDE  is  the  correlation  of  diagonal  errors  defined  by 
CDE=0.5*E{[X(p,q-l)-(XEST(p,q-l)]*[X(p-l,q)-XEST(p-l,q)] 

+  [  X  (  p  ,  q  )  -  XE  S  T  (  p  ,  q  )  ]  *  [  X  (  p  -  1  ,  q  -  1  )  -  XE  S  T  (  p  -  1  ,  q  -  1  )  ]  } 

(5.2) 

where  XEST  is  the  reconstructed  estimate  of  X,  and  E  is  the 
expectation  operator.  Here  E  may  be  thought  of  as  a  spatial 
operator  that  averages  over  all  pixel  locations  (p,q)  in  an 
image,  or  equivalently  it  may  be  thought  of  as  a  combination 
of  a  spatial  operator  averaging  over  an  N  by  N  block  of  an 
image  and  an  ensemble  operator  averaging  over  all  such 
blocks  in  an  image.  Both  viewpoints  will  be  used  in  the 
sequel . 

From  its  definition  the  CDE  is  seen  to  be  a  measure  of 
the  correlation  of  the  reconstruction  error  between  diagonal 
neighbors.  For  example,  if  a  hypothetical  transform 
resulted  in  a  reconstructed  image  that  was  merely  the 
original  image  decreased  by  5  intensity  units  at  all  pixel 
locations  (i.e.  XEST(p,q)  -  X(p,q)-5  for  all  (p,q)),  then 
diagonal  errors  would  be  totally  correlated  and  the  diagonal 
intensity  differences  upon  which  the  gradient  is  based  would 
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be  unaffected  (as  would  all  edge  information).  In  such  a 
hypothetical  case  IMSE*CDE»25  and  GMSE*0.  However,  as  seen 
in  Table  1,  the  Hotelling  transform  consistently  results  in 
negative  CDE  values,  causing  GMSE  to  be  greater  than  2*IMSE. 

Proof  of  Relationship: 

The  GMSE  is  the  expected  value  of  the  square  of  the 
norm  of  the  gradient  error  defined  by 

GMSE  -  E  { | |  G(p,q)  -  GEST(p,q)  )|**2}  (5.3) 
where  ||  ||  denotes  the  norm,  GEST(p,q)  represents  the 
estimate  of  the  gradient  at  location  (p,q)  obtained  from  a 
reconstructed  image,  and  the  expectation  operator  E  is 
defined  as  described  earlier.  This  may  be  written  making 
the  spatial  averaging  over  an  N  by  N  block  explicit  and 
reducing  E  to  a  mere  ensemble  expectation  as  shown  on  the 
next  page. 
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GMSE  = 


N  N 

£  £ 


4*N**2  p=l  q=l 
[  (X(p ,q-l)  -  XEST(p.q-l)) 


+  (X(p,q)  -  XEST(p.q)) 

-  (X(p-l,q-l)  -  XEST(p-l ,q-l) ) 

-  (X(p-l,q)  -  XEST(p-l,q))]**2 

+  [  ( X( p- 1 , q )  -  XEST(p-l.q)) 

+  (X(p,q)  -  XEST(p.q)) 

-  (X(p-l.q-l)  -  XEST(p-l.q-l)) 

-  (X(p ,q-l)  -  XEST(p,q-l))]**2 


N  N 

£  £ 

p=l  q=l 


4*N*  *2 


2 * [  (X(p,q-1)  -  XEST(p.q-l)  )**2 

+  (X(p,q)  -  XE ST (p,q))**2 
+  (X(p-l.q-l)  -  XEST(p-l,q-l))**2 
+  (X(p-l,q)  -  XEST(p-l,q))**2] 
-4*[  (X(p,q-1)  -  XEST(p,q-l)) 

*(X(p-l,q)  -  XEST(p-l.q)) 

+  (X(p,q)-XEST(p,q)) 
*(X(p-l,q-l)-XEST(p-l,q-l))]  J 

which  reduces  to 

GMSE  =  2* ( IMSE  -  CDE) 


(5.4) 


(5.5) 


Note  that  the  gradient  at  pixel  locations  in  the  top  row  and 
left  column  of  an  N  by  N  block  depend  on  neighboring  pixels 
above  and  to  the  left  of  the  block.  The  above  notation  is 
based  on  the  top  row  and  left  column  of  neighbors  being 


1  A 
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designated  row  0  and  column  0,  respectively. 

Within  minor  roundoff  errors  this  relationship  16 
verified  by  Table  1.  For  example  using  N-4  on  Image  A  gives 
CDE  =»  -76  and  IMSE  *  371  for  a  calculated  GMSE  of 
GMSE  -  2 * ( 3 7  1  +  76)  -  894 

which  is  reasonably  close  to  the  simulated  value  of  887. 


\ 


IS 


6.  A  GMSE  BASED  LINEAR  TRANSFORMATION 

The  gradients  of  an  N  by  N  block  of  pixel  locations 
within  an  Image  form  a  gradient  block  which  depends  on  the 


border 

pixels 

at  the 

top 

and  to  the  left 

of  the  N 

by  N 

block . 

A 

linear 

transformation  is 

introduced 

which 

individually 

encodes 

and 

data  compresses 

overlapping 

(N+l) 

by  (N+l)  intensity  blocks  so  as  to  minimize  the  resulting 
GMSE  for  the  N  by  N  gradients  which  are  calculated  for  each 
reconstructed  intensity  block.  This  formulation  allows  each 
N  by  N  gradient  block  to  be  independently  reconstructed 
using  the  transform  coefficients  generated  from  its 
corresponding  (N+l)  by  (N+l)  intensity  block.  As  with  the 


Hotelling 

transform,  the  (N+l) 

by 

(N+l) 

pixels  in 

each 

intensity 

block  are  row  scanned 

to 

form 

a  vector 

wi  th 

(N+l)**2  elements,  where  **  represents  exponentiation.  A 
complete  set  of  (N+l)**2  transform  coefficients  could  be 
obtained  by  multiplying  this  vector  by  an  (N+l)**2  by 
(N+l)**2  matrix  A.  To  achieve  data  compression  ISAVE 
coefficients  are  generated  (ISAVE  <  (N+l)**2)  using  an  ISAVE 
by  (N+l)**2  matrix  A.  To  minimize  IMSE  the  Hotelling 
transform  chooses  the  rows  of  A  to  be  eigenvectors  of  the 
image  block  covariance  matrix. 

To  determine  the  matrix  A  that  minimizes  GMSE  a  set  of 
summations  describing  the  GMSE  as  a  function  of  the  matrix  A 
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and  block  covariance  matrix  was  derived.  (This  equation  is 
derived  at  the  end  of  this  section  and  implemented  in  the 
attached  Fortran  77  Subroutine  DGMSE.)  For  a  given  block 
size,  block  covariance  matrix,  and  specified  data 
compression,  the  matrix  A  that  minimizes  GMSE  was  found 
numerically  using  a  steepest  descent  algorithm  based  on  an 
article  by  Fletcher  and  Powell  17].  This  algorithm  required 
the  derivation  of  an  equation  describing  the  gradient  of  the 
function  GMSE.  (This  equation  is  implemented  in  the 
attached  Subroutine  DGGMSE.) 

The  eigenvectors  used  by  the  Hotelling  transform  were 
used  to  form  an  initial  guess  for  the  matrix  A  that  would 
minimize  GMSE.  The  optimal  matrix  A  was  then  calculated 
using  the  Fletcher-Powell  search  technique.  For  several 
combinations  of  N  and  ISAVE,  Table  4  compares  the  resulting 
GMSE  using  the  Hotelling  eigenvector  matrix  to  the  GMSE 
achieved  using  the  optimal  matrix  A.  The  theoretical  values 
were  obtained  using  the  GMSE  equation  implemented  in 
Subroutine  DGMSE.  The  simulated  values  were  obtained  by 
actually  transforming,  data  compressing,  intensity 
reconstructing,  and  independently  calculating  the  gradients 
in  reconstructed  overlapping  (N+l)  by  (N+l)  intensity  blocks 
of  Image  A.  Compared  to  the  Hotelling  transform,  the 
optimal  matrix  A  resulted  in  18Z  to  52Z  lower  GMSE  in  the 
cases  studied. 

Calculation  of  the  optimal  matrix  A  using  the  iterative 
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FI e t che r-Po we  1 1  algorithm  is  quite  computationally 
intensive.  For  N=6  and  ISAVE»9,  a  moderate  block  size,  the 
matrix  A  consists  of  9*(6+l)**2  =  441  elements.  The 
Fletcher-Powell  algorithm  seeks  to  find  the  resulting 
441-dimensional  vector  that  minimizes  GMSE.  For  this 
particular  example  708  calls  to  the  search  routine  described 
in  [7]  were  needed,  consuming  approximately  200  hours  of  CPU 
time  on  a  Data  General  MV  10000  computer.  As  with  the 
Hotelling  transform  discussed  in  Section  4,  this 
disadvantage  may  be  overcome  by  storing  a  small  set  of 
selectable  pre-calculated  A  matrices  for  real-time  use. 

The  optimal  matrix  A  described  above  minimizes  the  GMSE 
of  gradient  blocks  in  an  image  that  are  reconstructed 
independently.  Similarly  an  optimal  matrix  A  calculated  for 
an  entire  image  (based  on  an  ensemble  covariance  matrix) 
without  partitioning  the  image  into  blocks  would  achieve 
minimum  GMSE  over  all  linear  transformations  (although 
calculating  A  for  entire  images  would  be  computationally 
prohibitive  using  iterative  search  techniques  such  as  the 
Fletcher-Powell  algorithm).  However  the  independent 
reconstruction  of  adjacent  gradient  blocks  does  not  allow 
correlations  between  adjacent  blocks  to  be  exploited.  Using 
the  Hotelling  transform  on  non-overlapping  N  by  N  intensity 
blocks  to  data  compress  and  reconstruct  an  entire  intensity 
image,  followed  by  gradient  calculations  on  the 
reconstructed  intensity  image  takes  advantage  of 
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block-to-block  correlations.  Table  5  compares  the  GMSE  for 
dependently  reconstructed  gradient  blocks  obtained  using  the 
Hotelling  transform  in  this  fashion  to  the  GMSE  obtained 
using  the  optimal  matrix  A  which  independently  reconstructs 
each  block.  For  the  combinations  of  N  and  ISAVE  shown  in 
Table  5,  adjacent  block  correlations  cause  the  GMSE  of 
independently  reconstructed  gradient  blocks  to  be  inferior. 
For  the  larger  block  sizes  shown  in  Table  5  the  advantage  of 
adjacent  block  correlations  seems  to  diminish.  For  block 
sizes  greater  than  N-6  it  is  not  yet  known  which  method  will 
result  in  lower  GMSE.  However,  as  stated  above,  if  the 
block  size  increases  to  encompass  the  entire  image,  the 
optimal  matrix  A  would  achieve  minimum  GMSE  while  the 
relative  performance  of  the  Hotelling  transform  is  unknown. 

Derivation  of  GMSE  as  a  function  of  A  and  CX: 

Equation  5.4  shows  that  the  GMSE  is  related  to  the 
expected  value  of  a  function  of  original  and  reconstructed 
block  pixel  intensity  values  X  and  XEST,  respectively.  At 
this  point  an  expression  is  derived  relating  GMSE  to  the 
elements  of  the  ISAVE  by  (N+l)**2  transform  matrix  A  and  the 
elements  of  the  (N+l)**2  by  (N+l)**2  block  covariance  matrix 
CX.  As  with  the  Hotelling  transform  in  Equation  4.1,  the 
ensemble  mean  for  each  block  pixel  is  subtracted  prior  to 
multiplication  by  the  transform  matrix  A.  Therefore  the 
following  derivation  can  be  simplified  by  assuming  the  image 
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block  values  X  have  zero  mean  which  reduces  Equations 
4. 1-4.3  to 


Y  =  AX 

(6.1) 

CX  =  E { XX ' } 

(6.2) 

XE  ST  -  A'Y 

(6.3) 

For  ease 

o  f 

presentation 

the  following 

derivation  is 

based  on  N=4 

so 

that  5  by 

5  overlapping 

blocks  are 

t  rans  f  ormed  . 

The 

extension  of 

the  results 

to  arbitrary 

values  of  N  is  trivial.  For  N=4,  Equation  6.1  may  be 
written  as 

N  N 

Y(i)  *£  £  A(i , 5m+n+l )*X(m,n)  (6.4) 

m=o  n=o 

where  Y(i)  is  the  ith  element  of  Y  for  i  from  1  to  23.  Then 
Equation  6.3  becomes 

I  SAVE 

XE  ST ( p  ,  q )  =  £  A(i,5p+p  +  l)*Y(i) 

i=l 

I  SAVE 

=  Y,  A(i,5p+q  +  l) 

i=l 

N  N 

*  £  £  A(i , 5m+n+l)*X(m,n)  (6.5) 

m=o  n=o 

for  p  and  q  from  0  through  4.  Expanding  the  quadratic  terms 
in  Equation  5.4  and  taking  the  expectation  of  each  resulting 
product  we  obtain  terms  of  the  following  forms: 
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E{X(p ,q)X(r ,8 ) }  -  CX( 5p+q+l , 5r+s+l ) 


(6.6) 


E(X(P.q)XEST(r,s)}  - 

I  SAVE 

-  H  A(i,5r+s+l) 
i-1 

N  N 

*1^0  rJo  A( i , 5m+n+l )  *  E { X(p , q ) X(m, n ) } 

ISAVE  N  N 

^  A(i , 5r+s+l)A(i , 5m+n+l) 
i*l  m*o  n*o 

*  CX( 5p+q+i , 5m+n+l )  (6  7 

E { XE  ST ( p , q ) XEST ( r , 8 ) }  - 
ISAVE  N  N 

E*  2-  2E  £  A(i , 5p+q+l)A(i , 5m+n+l)X(m,n) 

i*l  m=o  n=o 
ISAVE  N  N 

*  ^  2  E  A(i'  ,  5r+s  +  l)A(i '  ,5m'+n'+l) 

i’=l  m’*o  n'-o 

*  X(o ' ,  n '  )  } 

ISAVE  N  N  ISAVE  N  N 

-  IEEE  £  £  ( 

i“l  m”°  n-o  i'-l  m'-o  n’-o  1 

A(i,5p+q  +  l)A(i,  5m+n  +  l  )A(i  *  ,  5r+s+l  )A(  i  '  ,  5m '  +n  '  +1 ) 

*  CX( 5m+n+l ,  5m ' +n ' +1 )  I 

I  (6.8) 

sing  Equations  6.6  through  6.8  enables  us  to  now  rewrite 
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GMSE 


1 

4*N**2 


N  N 


E  £ 

p=i  q=i 


2[CX(5p+q+l,5p+q+l)+CX(5p+q,5p+q) 

+CX(5(p-l)+q  ,  5(p-l)+qHCX(5(p-l)+q  +  l,  5(p-l)+q  +  l)  ] 
-4[CX(5p+q, 5(p-l)+q+l)+CX(5p+q+l, 5(p-l)+q)) 


I SAVE  N  N 

•'III  A( i , 5m+n+l )  { 
i=l  m=o  n=*o 


A(i , 5p+q+l ) [CX( 5p+q+l , 5a+n+l) 

-CX( 5(p— 1) +q , 5m+n+l ) ] 

+A(i , 5p+q ) [CX( 5p+q , 5m+n+l ) 

-CX(5(p-l)+q+l, 5m+n+l ) ] 
+A(i,5(p-l)+q)[CX(5(p-l)+q,5n+n+l) 

-CX( 5p+q+l , 5m+n+l ) ] 

+A(i,5(p-l)+q+l) [CX(5(p-l)+q+l,5m+n+l) 

-CX( 5p+q , 5m+n+l ) ] 

I SAVE  N  N 

-0-5  £  E  Z 

i '  =1  m'=o  n'»o 

A(i* ,  5m'+n'  +  l)CX( 5m+n+l ,  5m'  +n ' +1 ) 

MA(i,5p+q  +  l)A(i'  ,5p+q  +  l)+A(i,5p+q)A(i'  ,5p+q) 
+A(i,5(p-l)+q)A(i' ,5(p-l)+q) 
+A(i,5(p-l)+q+l)A(i',5(p-l)+q+l) 
-2*A(i,5p+q)A(i' ,5(p-l)+q+l) 

-2*A(i , 5p+q+l )A(1 ' »5(p-l)+q)]  }(  (6.9) 
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i 


Equation  6.9,  derived  for  N  =  4,  can  be  generalized  by 
replacing  the  numeral  5  in  the  CX  and  A  matrix  indices  by 
(N+l).  This  generalized  expression  is  programmed  in  the 
attached  Subroutine  DGMSE  wherein  the  matrix  A  described 
here  is  row  scanned  to  form  a  vector  A  with  dimension 
ISAVE*(N+1)**2. 
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7. 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  Hotelling  transform  which  reduces  image  data  so  as 
to  minimize  intensity  mean-square  error  (IMSE)  in  the 
reconstructed  image  was  also  found  to  have  significantly 
better  edge  retaining  ability  than  simple  averaging.  The 
reconstructed  edges  were  quantitatively  compared  to  those  in 
the  original  images  using  the  MSE  based  and  receiver 
operating  characteristic  (PD  and  PF )  measures  described  in 
Section  3.  One  such  measure  used  was  the  gradient 
mean-square  error  (GMSE).  Both  the  reconstructed  IMSE  and 
GMSE  using  the  Hotelling  transform  tend  to  decrease  as  the 
encoding  block  size  increases.  An  equation  relating  GMSE  to 
IMSE  was  developed  in  Section  5. 

For  image  gradient  blocks  that  are  independently 
reconstructed.  Section  6  derives  the  linear  transformation 
matrix  A  that  minimizes  the  reconstructed  GMSE,  and  in  that 
sense  maximizes  edge  retention.  Calculation  of  the  optimal 
matrix  A  is  quite  computationally  intensive.  The  largest 
block  size  for  which  the  matrix  A  was  calculated  was  for 
overlapping  7  by  7  intensity  blocks  (N«6).  In  this  case  the 
GMSE  obtained  using  the  Hotelling  transform  on  overlapping  7 
by  7  intensity  blocks  to  independently  reconstruct  6  by  6 
gradient  blocks  of  Image  A  was  1042,  30Z  higher  than  the 
GMSE-601  obtained  using  the  optimal  matrix  A.  However, 
independent  reconstruction  of  gradient  blocks  does  not  allow 
correlations  between  adjacent  blocks  to  be  exploited. 
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Dependent  gradient  block  reconstruction  using  the  Hotelling 
transform  on  non-overlapping  6  by  6  intensity  blocks, 
reconstructing  the  resulting  intensity  image,  and  then 
applying  the  gradient  operator  resulted  in  GMSE=729,  9X 
lower  than  independent  reconstruction  using  the  optimal 
matrix  A.  Apparently  at  this  block  size  adjacent  block 
correlation  plays  a  dominant  role  in  determining  the 
reconstructed  GMSE. 

This  research  has  demonstrated  that  the  edge 
information  retained  in  a  data  compressed  image  can  best  be 
extracted  by  using  knowledge  of  the  image  data  reduction 
technique  used.  An  optimal  system  design  should  include 
selecting  the  image  reducing  technique  based  on  the 
reconstructed  end  product  desired.  For  example,  if  the 
desired  end  product  is  a  gradient  image,  the  intensity  image 
should  be  reduced  so  as  to  minimize  some  measure  of  the 
reconstructed  gradient  error  such  as  GMSE. 

Suggested  areas  for  future  research  include: 
Simulations  using  larger  block  sizes  comparing  independent 
and  dependent  gradient  block  reconstruction  errors. 

Development  of  a  transform  which  minimizes  GMSE  while 
including  the  effects  of  adjacent  block  correlation. 

Development  of  transforms  to  minimize  Klrsch  or  Sobel  edge 
errors. 

Development  of  transforms  for  maximum  edge  retention  based 
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on  specific  image  covariance  models. 

Determine  edge  retaining  abilities  of  other  non-image 
dependent  transforms. 
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Figure  1.  Test  Images  A, 


Table  1.  MSG  Based  Edge  Comparison  of  Averaging  vs 
Hotelling  Transform 


ROBERTS  KIRSCH 


IMSE 

GMSE 

MSE 

MSE 

CDE 

Image  A: 

Averaging : 

524 

1367 

1393 

520 

-164 

Hotelling 

N-4 

371 

887 

760 

366 

-76 

6 

305 

729 

641 

289 

-63 

8 

279 

704 

619 

261 

-78 

Image  B: 

Ave  raging : 

54 

146 

145 

36 

-19 

Hotelling 

N-4 

34 

84 

70 

29 

-8 

6 

26 

62 

54 

21 

-5 

8 

23 

54 

48 

18 

-4 

Image  C: 

Averaging : 

85 

226 

224 

59 

-29 

Hotelling 

N-4 

55 

134 

114 

50 

-13 

6 

43 

102 

92 

39 

-8 

8 

38 

91 

81 

34 

-8 
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Table  2 


Kirsch  Based  Edge  Comparison  of  Averaging  vs 
Hotelling  Transform 


THRESHOLDED 

KIRSCH 


MSE 

PD 

PF 

E 

Image  A: 

Original : 

0 

1 

0 

.780 

Averaging : 

231 

.469 

.0187 

.747 

Hotelling 

N»4: 

184 

.475 

.0071 

.791 

6: 

138 

.541 

.0061 

.784 

8  : 

116 

.569 

.0053 

.784 

PD  *  Probability  of  Detection 
PF  =  Probability  of  False  Alarm 
E  *  Average  Edge  Coherence 
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Table  3.  Sensitivity  of  GMSE  to  laage  Block 
Covariance  Matrix  CX 


N  -  4; 

I SAVE  -  4 

IMAGE 

IMAGE 

TRANSFORMED 

CX  USED 

GMSE 

A 

A 

887 

A 

B 

1090 

A 

C 

1098 

B 

A 

87 

B 

B 

84 

B 

C 

85 

C 

A 

139 

C 

B 

136 

C 

C 

134 

31 


Table  4.  GMSE  of  Independently  Reconstructed 
Gradient  Blocks  of  Image  A 


THEORETICAL 

SIMULATED 

GMSE 

GMSE 

N-2  ISAVE=2 : 

Hot  e 1 1 i ng : 

1391 

1374 

Optimal : 

914 

902 

N-3  ISAVE=4 : 

Hotelling : 

979 

965 

Optimal : 

727 

714 

N=4  ISAVE-4: 

Ho  tel ling : 

1242 

1227 

Optimal : 

1056 

1042 

N*6  ISAVE=9 : 

Hotelling : 

1042 

1030 

Optimal : 

801 

790 

32 


Table 


N 


N 


N 


N 


Independent  versus  Dependent  Gradient 
Block.  Reconstruction 


GMSE 

ISAVE=2 : 

Independent : 

902 

Dependent  : 

620 

ISAV£=4 : 

Independent : 

714 

Dependent  : 

424 

ISAVE=4 : 

Independent : 

1042 

Dependent : 

888 

I SAVE  =  9 : 

Independent : 

790 

Dependent  : 

729 
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APPENDIX 


1 


of 

Fortran  Subroutines 


SUBROUTINE  DGMSE(CX, A, N, ISAVE, IASIZE, ICXSIZ ,RMSE) 

C  CALCULATE  THEORETICAL  GRADIENT  MSE  TRANSFORMING  (N+1)X(N+1)  BLOCKS 
C  AND  SAVING  ISAVE  COEFFICIENTS. 

C  THE  VECTOR  A  IN  THIS  SUBROUTINE  IS  A  ROW  SCANNED  VERSION  OF 

C  THE  ORIGINAL  ISAVE  X  (N+l)  MATRIX  A. 

IMPLICIT  DOUBLE  PRECISION(A-H.O-Z) 

DOUBLE  PRECISION  CX( ICXSIZ , ICXSIZ) ,A( IASIZE) 

INTEGER  P,Q,U,V,C,D,E,F 

N1=N+1 

N1SQ=N1**2 

RMSE=O.DO 

DO  10  P=1,N 

DO  10  Q=1,N 

C=N1*(P-1)+Q 

D=C+1 

E=N1*P+Q 

F=E+1 

SUM1=2.D0*(CX(F,F)+CX(E,E)+CX(C, 

1C)+CX(D, D) )-4.D0* 

1 (CX(E,D)+CX(F ,C) ) 

SUM2=0.D0 

SUM3=O.DO 

DO  20  1=1, ISAVE 

U=( 1-1 )*N1SQ 

DO  20  M1=1,N+1 

DO  20  NN1=1 , N+l 

M=M1-1 

L=NN1-1 

MCX=N1*M+L+1 

SUM2=SUM2+ 

1A(U+F)*(CX(F,MCX)-CX(C,MCX))+ 

1A(U+E)*(CX(E,MCX)-CX(D,MCX))+ 

1A(U+C)*(CX(C,MCX)-CX(F,MCX))+ 

1A(U+D)*(CX(D,MCX)-CX(E,MCX) ) 

DO  30  J=l, ISAVE 
V=( J-l )*N1SQ 
DO  30  MP 1=1, N+l 
DO  30  NP1=1,N+1 
MP=MP1-1 
NP=NP1-1 
MPCX=N1*MP+NP+1 

SUM2=SUM2-0. 5D0*(A(U+F )*A( V+F )+ 

1 A( U+E ) *A( V+E ) +A( U+C ) *A( V+C )+ 

1A(U+D)*A(V+D)-2.D0* 

1A(U+E)*A( V+D)-2.D0*A(U+F)*A( V+C) )* 

1 A( V+MPCX)*CX(MPCX,MCX) 

30  CONTINUE 

SUM3=SUM3+SUM2*A(U+MCX) 

20  SUM2=0.D0 

10  RMSE=RMSE+SUM1-4.D0*SUM3 
RMSE=RMSE/(4.DO*DBLE(N**2) ) 

RETURN 

END 
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SUBROUTINE  DGGMSE(CX,A,G,N, ISAVE, IASIZE, ICXSIZ) 

C  CALCULATES  GRADIENT  OF  GMSE 

IMPLICIT  DOUBLE  PRECISION  A-H.O-Z) 

DOUBLE  PRECISION  G( IASIZE) , A( IASIZE) , CX( ICXSIZ , ICXSIZ) 

INTEGER  X,Y,P,Q,C,D,E,F,U,V 

N1=N+1 

N1SQ=N1**2 

X=1 

Y=0 

7  Y=Y+1 

IF (Y.GT.N1SQ)X=X+1 

IF(X.GT. ISAVE) RETURN 

IF(Y.GT.N1SQ)Y=1 

U=( X-I )*N1 SQ 

IC=0 

ID=0 

IE=0 

IF=0 

SUM=0.D0 
DO  10  P=1 , N 
DO  10  Q=1,N 
C=N1*(P-1 )+Q 
D=C+1 
E=N1*P+Q 
F=E+1 

SUMJMN=O.DO 
DO  15  J=l, ISAVE 
V=( J-l )*N1SQ 
DO  15  MP=0,N 
DO  15  NP=0, N 
MPCX=N1*MP+NP+1 

1 5  SUMJMN=SUMJMN+(A(U+C)*A(V+C)+A(U+D)*A(V+D)+A(U+E)*A(V+E)+ 
1A(U+F)*A(V+F)-2.D0*(A(U+E)*A(V+D)+A(U+F)*A(V+C) ) )* 

1A( V+MPCX)*CX(MPCX, Y) 

10  SUM=SUM-0.5D0*SUMJMN+A(U+C)*(CX(C,Y)-CX(F,Y))+ 

1A(U+D)*(CX(D,Y)-CX(E,Y))+A(U+E)*(CX(E,Y)-CX(D,Y))+ 
1A(U+F )*(CX(F , Y)-CX(C,Y) ) 

SUMMN=0.D0 

DO  20  M=0,N 

DO  20  L=0,N 

SUMCX=0.D0 

JCX=N1*M+L+1 

ICNT=0 

JCNT-0 

DO  25  JY-l.N 
26  ICNT-ICNT+1 
JCNT-JCNT+1 
IF( ICNT.GT.N)GO  TO  25 
IF( Y.EQ. JCNT)GO  TO  30 
GO  TO  26 


25  ICNT=0 
GO  TO  35 
30  IC=1 

SUMCX=SUMCX+CX( Y, JCX)-CX( Y+Nl+1 , JCX) 
35  ICNT=N1+1 
JCNT=N*N1+1 
DO  45  JY=1 ,N 
46  ICNT=ICNT-1 
JCNT=JCNT-1 
IF( ICNT.EQ. 1 )G0  TO  45 
IF( Y.EQ. JCNT)GO  TO  50 
GO  TO  46 
45  ICNT=N1+1 
GO  TO  55 
50  ID=1 

SUMCX=SUMCX+CX( Y , JCX) -CX( Y+N 1 - 1 , JCX) 
55  ICNT=0 
JCNT=N1 
DO  65  JY=1,N 
66  ICNT=ICNT+1 
JCNT=JCNT+1 
IF(ICNT.GT.N)GO  TO  65 
IF(Y.£Q.JCNT)GO  TO  70 
GO  TO  66 
65  ICNT=0 
GO  TO  75 
70  IE=1 

sumcx=sumcx+cx(y,jcx)-cx(y-ni+i,jcx) 

75  ICNT=N1+1 
JCNT=N1SQ+1 
DO  85  JY=1 ,N 
86  ICNT=ICNT-1 
JCNT=JCNT-1 
IF ( ICNT.EQ. 1 )GO  TO  85 
IF( Y.EQ. JCNT)GO  TO  90 
GO  TO  86 
85  ICNT=N1+1 
GO  TO  20 
90  IF=1 

SUMCX=SUMCX+CX(Y, JCX)-CX( Y-Nl-1 , JCX) 
20  SUMMN=SUMMN+A( U+JCX) *SUMCX 
SUM=SUM+SUMMN 
SUMIMN=0.D0 
DO  100  1=1 , ISAVE 
V=( 1-1 )*N1SQ 
DO  100  M=0,N 
DO  100  L=0,N 
SUMPQ=0.D0 
JCX=N1*M+L+1 
DO  110  P=1 ,N 
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do  no  q=i,n 

C=N1*(P-1)+Q 

D=C+1 

E=N1*P+Q 

F=E+1 

110  SUMPQ=SUMPQ+CX(Y, JCX)*( 

1A(V+C)*A(U+C)+A(V+D)*A(U+D)+A(V+E)*A(U+E)+A(V+F)*A(U+F)-2.D0* 
1 (A(V+E)*A(U+D)+A(V+F)*A(U+C) ) ) 

100  SUMIMN=SUMIMN+A(V+JCX)*SUMPQ 
S(JM=SUM-0.  5D0*SUMIMN 
SUMMNP=0.D0 
ACNT=DBLE(IC+ID+IE+IF) 

DO  200  M=0,N 
DO  200  L=0,N 
DO  200  MP=0,N 
DO  200  NP=0,N 
SUMA=0.D0 
MCX=N 1 *M+L+1 
MPCX=N1*MP+NP+1 
SUMI=0.D0 
DO  210  1=1 , ISAVE 

210  SUMI=SUMI+A((I-1)*N1SQ4-Y)*A((I-1)*N1SQ+WCX) 
SUMA=SUMA4A(U+MPCX)*SUMI 
SUMJ=0.D0 
DO  220  J=l, ISAVE 

220  SUMJ=SUMJ+A( (  J-1)*N1SQ+-Y)*A(  (  J-l  )*N1SQ+MPCX) 

S UMA=SUMA+A( U4MCX) *SUMJ 
S  UMA=S  UMA*ACNT 
IF(IC.EQ.0)G0  TO  259 
SUMI=0.D0 
DO  250  1=1, ISAVE 

250  SUMI=SUMI+A( (I-1)*N1SQ+Y+N1+1 )*A( (1-1 )*N1SQ+MCX) 
SUMA=SUMA-2.D0*A(U+MPCX)*SUMI 

259  IF(ID.EQ.O)GO  TO  269 
SUMI=O.DO 

DO  260  1=1, ISAVE 

260  SUMI=SUMI+A( ( 1-1 )*N1S<*FY+N1-1)*A( (1-1 )*N1SQ+«CX) 

SUMA=SUMA-2 . DO*A(U+MPCX)*SUMI 

269  IF(IE.EQ.O)GO  TO  279 
SUMJ=0.D0 

DO  270  J=l, ISAVE 

270  SUMJ=SUMJ+A( ( J-l )*NlSQfY-Nl+l )*A( (J-l )*N1SQ+MPCX) 
SUMA=SUMA-2.D0*A(U+MCX)*SUMJ 

279  IF( IF .EQ.O)GO  TO  200 
SUMJ-O.DO 

DO  280  J-l, ISAVE 

280  SUM J-SUM J+A(  (  J- 1 )  *N  1  SQfY-N  1 - 1 ) *A( ( J- 1 ) *N 1  SCfr+MPCX) 

SUMA-SUMA-2 . DO*A( U4MCX)*SUMJ 

200  SUMMNP-SUMMNP+SUMA*CX(MPCX,MCX) 

SUM-SUM-0. 5D0*SUMMNP 
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